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1.  Introduction 

The  numerical  solution  of  a  wave  propagation  problem  in  a  very  large  or  unbounded  domain  provides  a  challenging  com¬ 
putational  difficulty  -  namely,  solving  the  problem  on  a  finite  computational  domain  while  maintaining  the  true  essence  of 
the  solution.  One  of  the  modern  techniques  that  has  garnered  a  significant  amount  of  attention  in  handling  this  challenge  is 
the  absorbing  or  non-reflecting  boundary  condition  (NRBC)  method.  In  using  this  method,  the  original  infinite  domain  is 
truncated  by  an  artificial  boundary  5,  resulting  in  a  finite  computational  domain  Q  and  the  residual  domain  D. 

When  truncating  the  domain,  the  modeler  must  devise  boundary  conditions  for  the  truncated  domain.  Of  course,  by 
imposing  a  boundary  where  one  does  not  physically  exist,  the  problem  is  changed  -  and  unless  chosen  carefully,  would  cer¬ 
tainly  be  expected  to  pollute  the  solution  as  the  problem  evolves  and  impinges  on  the  boundary.  Therefore,  two  main  pos¬ 
sibilities  exist  for  the  modeler: 

•  Choose  a  convenient,  easily  implementable  boundary  condition  that  does  not  necessarily  reflect  the  physical  problem  and 
solve  it  on  a  large  sub-domain.  The  idea  behind  this  technique  is  that  the  boundary  effects  are  negligible  for  a  short  time 
evolution  of  the  problem  in  a  small  area  of  interest  away  from  the  boundaries. 

•  Choose  a  boundary  condition  that  preserves  the  true  behavior  of  the  infinite  solution  at  the  boundary  and  solve  the  prob¬ 
lem  on  a  smaller  sub-domain.  The  idea  behind  this  technique  is  that  the  additional  effort  extended  to  impose  a  better 
boundary  condition  will  be  worth  the  effort  and  allow  for  solving  the  problem  on  a  smaller  domain. 

For  obvious  reasons,  the  first  possibility  has  only  limited  usefulness.  To  see  why,  suppose  that  we  wanted  to  model  the 
wave  motion  following  a  pebble  dropped  in  the  center  of  a  large,  still  pond.  Now,  suppose  that  we  have  a  truncated  domain 
to  model  this  phenomena  -  say  a  bathtub.  If  the  pebble  is  dropped  in  the  bathtub,  the  waves  generated  by  the  pebble  would 
propagate  much  like  that  in  the  pond  -  until,  that  is  -  the  wave  front  reaches  the  hard  walls  of  the  bathtub.  At  this  point,  the 
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bathtub  model  ceases  to  be  a  useful  representation  of  the  pond  due  to  behavior  caused  by  the  non-physical  boundary.  If  the 
modeler  wishes  to  see  what  happens  a  short  time  later  -  a  larger  bathtub  would  be  required.  This  same  principle  would  ap¬ 
ply  for  the  numerical  solution  of  this  propagation  problem  -  a  poor  choice  of  boundary  condition  mandates  the  use  of  a  lar¬ 
ger  computational  domain.  This  in  turn,  requires  additional  computational  resources.  For  this  reason,  much  effort  has  and 
continues  to  be  exerted  on  finding  stable,  efficient,  accurate  and  practical  means  of  reducing  this  reflection  through  so-called 
NRBCs  [1]. 

Several  high-order  NRBCs  have  been  devised  to  reduce  spurious  reflections  that  would  pollute  the  solution.  Beginning  in 
the  late  1980’s,  the  well-known  Engquist-Majda  [2]  and  Bayliss-Turkel  conditions  [3]  gave  way  to  Collino’s  [4]  low  deriv¬ 
ative,  auxiliary  variable  formulation  for  the  2D  scalar  wave  equation.  This  sparked  a  flurry  of  activity  in  an  effort  to  And  qual¬ 
ity,  high-order  NRBCs  that  were  easily  implementable.  See  [5-8]  for  reviews  on  the  subject. 

The  starting  point  for  the  family  of  NRBCs  discussed  here  is  the  condition  devised  by  Higdon  in  a  series  of  papers  [9-14], 
which  was  demonstrated  in  a  low-order  finite  difference  setting.  While  in  theory,  Higdon’s  NRBC  is  considered  a  high-order 
NRBC,  the  formulation  requires  evaluation  of  increasing  high-order  spatial  and  temporal  derivatives  as  the  order  of  the  NRBC 
increased.  In  [15]  Givoli  and  Neta  directly  extended  the  Higdon  scheme  to  high-order  finite  difference  discretizations  via  an 
algorithm  where  the  order  of  the  NRBC  was  simply  an  input  parameter.  Long  term  stability  using  this  formulation  was  dem¬ 
onstrated  in  [16].  They  later  extended  this  formulation  to  one  that  does  not  involve  any  high  derivatives  (hereafter  referred 
to  as  the  G-N  formulation).  The  elimination  of  all  high-order  derivatives  is  enabled  through  the  introduction  of  special  aux¬ 
iliary  variables  on  B.  This  construction  demonstrated  in  [17,18]  for  finite  differences  was  further  extended  in  [19]  for  finite 
element  schemes  to  solve  the  dispersive  wave  equation. 

Hagstrom  and  Warburton  [20]  (H-W)  also  used  the  Higdon  and  auxiliary  variable  framework  to  develop  a  symmetric 
boundary  formulation  in  a  full-space  configuration  where  special  corner  compatibility  conditions  were  developed  for  the 
non-dispersive  wave  equation.  One  reviewer  pointed  out  a  recently  published  article  by  Becache  et  al.  [21]  where  the  H- 
W  formulation  is  extended  for  anisotropic  and  convective  wave  equations  in  a  finite  element  setting.  Results  are  presented 
for  the  anisotropic  case  demonstrating  significant  improvements  as  the  order  of  the  NRBC  is  increased,  although  convergence 
appears  to  be  limited  by  the  numerical  discretization. 

Several  extensive  comparisons  of  the  G-N  and  H-W  formulations  [21-23]  have  shown  certain  advantages  of  the  H-W  for¬ 
mulation.  One  such  advantage  of  the  H-W  formulation  is  that  the  auxiliary  variables  are  shown  to  decay  as  the  order  of  the  NRBC 
is  increased.  While  this  property  is  not  required  for  convergence,  it  is  a  useful  property  if  dynamic  adjustment  of  the  NRBC  order 
is  desired.  Additionally,  while  both  the  G-N  and  H-W  formulations  decay  exponentially  fast  as  a  function  of  the  NRBC  order,  the 
H-W  formulation  has  a  reflection  coefficient  that  decays  twice  as  fast  as  the  G-N  formulation  [22,  p.  3676].  On  the  other  hand, 
the  G-N  formulation  is  convenient  in  that  through  a  judicious  choice  of  NRBC  parameters,  the  G-N  formulation  can  be  reduced 
to  a  first-order  (in  time)  system  which  is  extremely  sparse  in  nature  and  quite  easy  to  implement.  This  choice  of  NRBC  param¬ 
eters  does  not  hamper  convergence  properties  of  the  boundary  scheme. 

To  see  how  these  auxiliary  variable  formulations  are  constructed,  consider  Fig.  1  that  illustrates  the  NRBC  set-up  using  an 
infinite  wave-guide.  Here,  the  artificial  boundary  B  extends  from  the  southern  (Ts)  to  the  northern  {r^)  boundaries  of  the 
wave-guide,  thus  creating  the  east  (Fe)  and  west  (Fw)  boundaries  of  ^2  at  x  =  Xw  respectively.  Outside  of  the  area  enclosed 
by  these  boundaries  is  the  residual  infinite  domain  D. 

Once  a  suitable  NRBC  is  devised,  the  problem  is  solved  numerically  in  ^2,  by  the  finite  difference,  or,  as  in  the  case  of  this 
analysis,  the  spectral  element  (SE)  method.  The  SE  method,  originally  introduced  by  Patera,  . .  combines  the  generality  of 
the  finite  element  method  with  the  accuracy  of  spectral  techniques. . [24,  p.  468].  This  generalized  high-order  finite  ele¬ 
ment  method  selects  the  integration  and  interpolation  points  carefully  in  order  to  yield  accurate  but  efficient  solutions. 
See  e.g.  Giraldo-Restelli  [25]  for  more  details  on  this  method. 

In  the  present  paper,  the  G-N  auxiliary  formulation  of  the  Higdon  NRBC  is  implemented  on  a  reduced  form  of  the  line¬ 
arized  shallow  water  equations  (SWE)  under  non-zero  advection.  This  formulation  has  been  previously  demonstrated  in  a 
finite  difference  formulation  to  arbitrarily  high  NRBC  order  [26],  however,  accuracy  gains  realized  by  increasing  the  NRBC 
ceased  after  order  2.  The  formulation  used  here  seeks  to  remedy  this  limitation  by  using  a  high-order  treatment  of  space 
(SE)  and  time  (Runge-Kutta)  to  show  the  benefits  of  using  the  high-order  boundary  (G-N)  scheme.  Specifically,  the  interior 
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Fig.  1.  An  infinite  wave-guide  truncated  by  artificial  boundaries  Fw  and  Fe. 
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and  boundary  formulations  are  discretized  using  high-order  basis  functions  in  a  stable,  equal-order  interpolation  scheme  for 
all  the  variables  (this  does  not  violate  the  inf-sup  condition).  High-order  time  integration  is  performed  as  well  in  an  effort  to 
balance  all  of  the  errors  involved  with  the  numerical  solution.  The  computational  effort  associated  with  the  high-order 
boundary  scheme  can  be  shown  to  grow  only  linearly  with  the  order  [19]. 

It  should  be  noted  that  the  only  other  spectral  element,  high-order  boundary  approach  (to  the  authors’  knowledge)  are 
Kucherov  and  Givoli  [27]  and  Lindquist  et  al.  [28].  Kucherov  and  Givoli  demonstrate  exponential  error  convergence  of  the 
classical  wave  equation  on  a  semi-infinite  channel  when  solved  using  spectral  elements  and  high-order  boundary  treatment 
(using  the  H-W  boundary  scheme).  They  further  show  how  the  spectral  element  formulation  allows  the  NRBC  to  realize  its 
true  potential,  prior  masked  by  low-order  numerical  schemes.  They  note  that,  “Although  it  is  generally  felt  that  there  is  no 
need  to  treat  the  time  domain  ‘spectrally’  like  the  spatial  domain,  [they]  feel  that  a  consistently  high-order  treatment  re¬ 
quires  that  the  entire  approximation  be  spectral,  i.e.  the  convergence  of  all  three  types  of  error  -  the  spatial  and  temporal 
discretization  errors  and  the  ABC  error  -  be  exponential.”  [27] 

Subsequent  spectral  element  and  boundary  work  by  the  current  authors  [28]  using  the  G-N  formulation  of  the  dispersive 
wave  equation  on  a  semi-infinite  channel  showed  similar  results  to  those  presented  by  Kucherov  and  Givoli  and  how  a  high- 
order  treatment  of  the  time  domain  (up  to  order  10)  produces  additional  improvements.  These  improvements,  however, 
have  their  limits  thus  confirming  the  hypothesis  of  Kucherov  and  Givoli.  The  key  difference  in  this  work  is  that  we  extend 
the  high-order  space,  boundary  and  time  integration  results  previously  demonstrated  in  a  non-zero  advection  setting  to  one 
where  the  wave  medium  is  not  at  rest. 

The  following  is  the  outline  of  the  rest  of  this  paper.  In  Section  2,  the  problem  under  investigation  is  stated.  In  Section  3, 
an  overview  of  the  G-N  auxiliary  formulation  is  presented.  In  Section  4,  a  SE  semi-discrete  formulation  that  incorporates  the 
G-N  NRBC  with  any  desired  order  is  constructed.  In  Section  5  the  time-integrator  used  to  March  the  equations  in  time  is 
discussed.  The  performance  of  the  method  is  demonstrated  in  Section  6  via  numerical  examples. 

2.  Statement  of  the  problem 

To  motivate  the  problem  under  consideration,  consider  the  SWE: 

dtU  +  UdxU  +  vdyU  -fv=  -gdxf], 

dtV^udxV+  vdyV+fu  =  -gdyf],  (1) 

dtt]  +  dx(Hu)  dy(Hv)  =  0. 

We  use  the  following  shorthand  for  partial  derivatives 


The  shallow  water  model  in  its  current  form  is  non-linear.  We  have  three  state  variables:  u{x,y,t)  and  v{x,y,t)  are  the  un¬ 
known  velocities  in  the  x  and  y  directions  and  f/(x,y,t)  is  the  water  depth  above  a  reference  value.  Eurther,  H  is  the  water 
depth  as  shown  in  Eig.  2  such  that  H  =  hg  +  ^,/is  the  Coriolis  parameter,  and  g  is  the  gravity  acceleration.  Now,  suppose  that 
the  bottom  topography  is  flat  such  that  Hb  is  constant  and  u  and  z/can  be  described  by  a  constant  mean  term  and  a  small  0{S) 
deviation  from  that  value,  i.e. 

u  =  U-\-u*  v  =  V+v*. 

To  be  clear,  U  and  V  are  the  mean  velocities  with  respect  to  the  coordinate  axes.  Using  these  substitutions  and  neglecting  any 
0(3^)  terms  results  in  the  linearized  form  of  the  SWE: 
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dtU*  +  VdxUr  +  VdyU*  -f(V  +  P*)  =  -gdxf], 

dtP*  +  UdxP*  +  VdyP*  +  f  {U  +  u*)  =  —gdyf],  (2) 

dtL]  +  UdxL]  +  V  dyf]  +  +  dyP*^  =  0. 

van  Joolen  shows  in  [30],  how,  using  the  operator: 

—  =  —  +  U—  +  V—, 

Dt  dt  dx  dy’ 

the  SWE  can  be  reduced  to  a  single  variable  Klein  Gordon  equation  (KGE)  equivalent  for  the  wave  perturbation  vf  under  non¬ 
zero,  constant  advection  velocities  U  and  V.  Once  f]  has  been  found,  the  unknown  velocities  can  be  computed  in  a  similar 
manner  (see  [29]  for  details).  This  work  seeks  to  numerically  solve  this  two-dimensional  advective,  dispersive  wave 
equation: 

d^r,  +  (U^  -  Cl)d^,r,  +  (V^  -  Cl)dlr,  +  2Udl^  +  2Vdl^  +  2{Nd%^  +fr,  =  0,  (3) 

where  Cq  =  gH,  using  continuous  Galerkin  methods. 

While  the  assumptions  for  this  reduced  form  of  the  SWE  may  undermine  the  predictive  capability  of  the  already  simpli¬ 
fied  equations  of  fluid  motion,  the  advective  KGE  serves  as  an  important  test-bed  to  extend  and  improve  non-refiecting 
boundary  conditions  for  wave  propagation  in  a  non-static  environment.  When  shown  to  improve  performance  in  test  cases 
as  those  presented  here,  these  conditions  will  then  be  extended  to  the  full  linearized  SWE  system  and  further  to  include  non¬ 
linear  effects. 

We  begin  this  analysis  by  considering  the  semi-infinite  wave-guide  where  Dirichlet  boundary  conditions  17(0,3/,  t)  =f(y,t) 
are  prescribed  on  Tw,  no-fiux  boundary  conditions  0  =  0  on  Tjv  and  Fs  and  the  G-N  NRBCs  are  enforced  on  Fe.  Eurther,  ini¬ 
tial  data  rj{x,y,0)  and  dtrj{x,y,0)  are  prescribed  at  t=  0.  The  set-up  is  illustrated  in  Eig.  3. 


3.  G-N  auxiliary  variable  formulation 


We  present  a  brief  summary  of  the  G-N  auxiliary  variable  process  as  described  in  [19].  This  auxiliary  formulation  begins 
with  the  Higdon  [14]  boundary  condition  given  by: 


Hj: 


r]  =  0  on  Te. 


(4) 


Note  that  the  first-order  condition  7=1  corresponds  to  the  Sommerfeld  radiation  condition,  provided  that  Q  =  Co.  Auxiliary 
functions  (/)i,. .  .,7)/_i,  which  are  defined  on  Fe  as  well  as  in  the  exterior  domain  D  (see  Eig.  3)  are  now  introduced.  Eventually 
we  shall  use  these  functions  only  on  Fe,  but  the  derivation  requires  that  they  be  defined  in  D  as  well,  or  at  least  in  a  non¬ 
vanishing  region  adjacent  to  Fe.  The  functions  4>j  are  defined  via  the  relations 


II 

(5) 

(N 

II 

(6) 

+  =0.  (7) 

By  definition,  these  relations  hold  in  D,  and  also  on  Fe.  It  is  easy  to  see  that  (5)-(7),  when  imposed  as  boundary  conditions  on 
Fe,  are  equivalent  to  the  single  boundary  condition  (4).  If  we  also  define 


Fn 


Fig.  3.  A  semi-infinite  channel  domain  Q  truncated  by  artificial  boundary  Fe. 
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5 


(8) 


then  we  can  write  (5)-(7)  concisely  as: 

+  =./.j  j  =  (9) 

This  set  of  conditions  involves  only  first-order  derivatives.  However,  due  to  the  appearance  of  the  x-derivative  in  (9),  one 
cannot  discretize  the  (/>/  on  the  boundary  Fe  alone.  Therefore,  we  shall  manipulate  (9)  in  order  to  get  rid  of  the  x-derivative. 

The  function  rj  satisfies  the  dispersive,  advective  wave  Eq.  (3)  in  D.  Since  the  function  is  obtained  by  applying  the  linear 
(constant  coefficient)  operator  (^dx  +  to  it  is  clear  that  4>^  should  also  satisfy  the  same  equation  in  D.  Further,  since  (/)y 
is  obtained  by  applying  the  same  linear  operator  j  -  1  times  to  (/>!,  the  functions  4>j  should  satisfy  an  equation  like  Eq.  (3), 
namely, 

9^^.  +  (u^  _  +  (v^  -  +  2ml4>i  +  2Vdl4>i  +  2UVd%^j  +f4>i  =  0.  (10) 

Using  the  following  identities: 

dl^j  =  (^5.  -  ^5,)  (^9.  +  ^ ipj, 

dlt4>i  =  dt{dx4>i), 
dly<Pj  =  dy{dx<pj) 

and  combining  with  (10)  allows  us  to  write  (9)  as: 


+  (^(u^  -  C^)  (^  +  ^)  -  =  (u^  -  for  J  =  1, ...  J  -  1.  (11) 


In  (11)  and  elsewhere,  a  prime  indicates  differentiation  with  respect  to  y  along  r,  i.e.  the  tangential  derivative  along  F.  As 
desired,  the  new  boundary  condition  (11)  does  not  involve  x-derivatives.  In  addition,  there  are  no  high-y  or  t  derivatives  be¬ 
yond  second  order.  It  should  be  noted  that  van  Joolen  et  al.  [26]  developed  an  equivalent  formulation  using  Lagrangian 
derivatives. 

Rewriting  (5),  (11)  and  (8),  the  new  formulation  of  the  Jth-order  NRBC  on  F  can  be  summarized  as  follows: 


M  +  dxn  =  h, 

(12) 

+ 

II 

1 

1 

-t- 

1 

+ 

(13) 

o 

III 

III 

O 

(14) 

Po  =  f,  =  Cj  =  ^-2V,  Xy  =  V^-Cl 

Li  Lj  Cj  Cj 

p.=  (u^-Cl)Q-  +  ^'^-2U,  y  =  2UV,  A.  =  l/'-C^. 

4.  Spectral  element  method 

The  spectral  element  (SE)  method  is  a  generalized  high-order  finite  element  method  where  the  integration  and  interpo¬ 
lation  points  are  selected  carefully  in  order  to  yield  accurate  but  efficient  solutions.  For  this  problem,  we  will  discuss  two 
formulations-one  for  the  interior  and  one  for  implementing  the  boundary  conditions. 

4. 1 .  Interior  formulation 

The  weak  form  of  the  problem  in  Q  is  constructed.  The  solution  is  sought  in  the  space  of  trial  functions, 

V  =  {rj\rj  G  H^(0)  and  rj  =  0  on  Fw}.  (15) 

Here,  H\^2)  is  the  Sobolev  space  of  functions  whose  first  derivatives  are  also  square-integrable  in  Q.  Now,  Eq.  (3)  is  multi¬ 
plied  by  the  basis  functions  *F,(x,y)  g  V  and  integrated  over  Q  so  the  weak  form  is: 
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[  +  4  [  Widlt]dQ  +  2.y  [  Widlf]dQ  +  2U  [  Wid^i]dQ  +  2V  [  WidyijdQ  +  UV  (  Wid^  t]dQ 

Jq  Jq  Jq  Jq  Jq  Jq 

+  UV  [  Wid^  t]dQ+p  [  >Pit]dQ  =  0. 

Jq  Jq 


(16) 


In  order  to  maintain  a  symmetric  form  of  the  problem,  the  mixed  derivative  has  been  split  appropriately.  To  ensure  the 
solution  vf  is  in  the  vector  space  requires  some  special  handling  of  the  second  order  derivatives,  which  is  facilitated  by  the 
use  of  Green’s  theorem,  i.e. 


/  <Pid^^t]dQ=  [  <Pid,t]n,dr-  [  d,Wid,t]dQ, 

J  Q  Jr  J  Q 

where  n  is  the  outward  normal  on  r  and  is  the  x-component  of  that  outward  normal. 

Extending  this  idea  for  each  second  order  derivative  in  (16)  gives  the  weak  form: 

[  Wii^dQ-1^  [  d^Wid^rjdQ-^y  [  dyWidyfjdQ  -  UV  [  dyWid^rjdQ  -  UV  [  d^WidyfjdQ +  2U  [  Wid^fjdQ 
Jq  Jq  Jq  Jq  Jq  Jq 

+  2V  (  Widyi^dQ+f  [  Wit]dQ  +  X^  (  Wid^rjn^dr  +  Xy  (  Widyt]nydr +  UV  [  Wid.rjnydr 
Jq  Jq  Jr  Jr  Jr 


(17) 


UV 


Widytjn^ 


dr  =  o. 


(18) 


Using  the  no-flux  boundary  conditions  on  the  northern  and  southern  borders  along  with  the  normal  vectors  specified  by 
the  structured,  rectangular  grid  shown  in  Fig.  4,  the  problem  may  be  simplified  to  account  for  contributions  on  individual 
boundaries.  Using  this  information  along  with  (12)  to  eliminate  the  normal  derivative  term  on  r^,  we  get: 

[  Wii^dQ-2^  [  d^Wid^rjdQ-2y  [  dyWidyfjdQ  -  UV  [  dyWid^rjdQ  -  UV  [  d.WidyrjdQ +  2U  [  Wid^fjdQ 
Jq  Jq  Jq  Jq  Jq  Jq 

+  2V  f  Widyt]dQ+p  [  Wit]dQ  +  X,  [  Wi(l>,drE-lioxU  WitjdrE  +  UV  [  Wid,t]drN-UV  [  Wid,t]drs 

J  Q  J  Q  J  Tp  J  r  p  J  r  fj  J  r  Q 


UV  [  WidyrjdrE  =  0. 
Jle 


(19) 


4.2.  Boundary  formulation 


Since  the  term  appears  in  the  interior  formulation  (19),  this  is  not  a  complete  weak  form  of  the  problem  in  Q.  We  turn 
our  attention  to  computing  a  separate  weak  form  for  (13)  to  complete  the  problem  statement. 

As  in  the  interior  formulation,  we  multiply  (13)  by  the  test  function  tj  and  integrate  it  over  Fe.  After  integration  by  parts 
and  simplifying  yields: 

(Xj  I  dr£  +  {j  j  Xi^p^dre  +  Xy  I  drE  + /Jj  j  x^jdCE-y  I  Xip^drE-p  I  Ti</>j_idrE 

dTp  dr^  dTp  dr^  dr^  dr^ 

=  4  /  dFE  (20) 

Jte 

for  j  =  1,. .  .,7  -  1,  e  H^iFE)  and  any  t/  e  H^iFE).  Recall  from  the  auxiliary  variable  formulation  (14)  that  (/)o  =  ^  and  =  0. 


n  = 


O' 

1 
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The  formal  problem  statement  is  then:  Find  rj  eV  and  (/)j  g  =  1,. . .  J  -  1,  such  that  Eqs.  (19)  and  (20)  are  satisfied 

G  V  and  t,  g  H\rE). 

43.  Galerkin  expansion 


We  now  turn  our  attention  to  the  spatial  discretization.  First,  we  expand  the  solution  variables  rj  and  (/)y  using  the  same 
basis  functions  used  in  the  weak  form  as  follows: 


=  =  j=  1,2,... J-1. 


(21) 


Here,  Np  refers  to  the  number  of  points  that  Q  is  discretized  into  and  refers  to  the  number  of  points  that  Fe  is  discretized 
into.  Next,  we  substitute  this  basis  function  expansion  directly  into  the  weak  form,  resulting  in  the  following  matrix  form  of 
this  problem: 


Mp  +  (2UDx  +  2VDy^i]  +  {~2xLxx  —  2yLyy  —  UVLyx  —  UVLxy^f]  +  2xAe(I)\  —  (^o^xEeLJ 

+  {UVCm  -  UVCs  +  UVCe)!]  =  0, 


(22) 


+  (XyL”  -fM”)  =  X,M'’<Pjy,  j  =  1 1 ,  (23) 

where  M,  Dx,  Dy,  L  xx^  Lyy,  Lxy  and  Lyx  are  interior  formulation  matrices  of  size  Np  x  Np.  Further,  Be,  Cjv,  C5  and  Ce  are  interior 
formulation  matrices  integrated  along  the  boundaries  of  size  Np  x  Np  while  Ae  is  of  size  Np  x  N^.  Finally,  M^,  and  are 
auxiliary  formulation  matrices  of  size  Nt,  x  Nt,.  These  global  matrices  are  obtained  from  analogous  element  arrays  via  assem¬ 
bly,  given  by: 


Ne  Ne  Ne  Ng 

M=/\M^  Be  =  /\K,  D,  =  /\DI  Dy  =  /\D;, 

e=l  e=l  e=l  e=l 


Ne 

=  AC 

^yy 

II 

Ne 

c  =  AC 

-3 

II 

e=l 

Nb 

e=l 

Nb 

e=l 

Nb 

e=l 

Ne 

(24) 

M”  = 

b=\ 

Ne 

=  Ad^ 

b=\ 

Ne 

L^  =  A'^ 

b=\ 

Ne 

<ii 

II 

Cs  =  AC 

Ae  =  /\AI, 

e=l  e=l  e=l 


where  is  the  direct  stiffness  summation  operator  required  by  all  continuous  Galerkin  methods.  The  expressions  for  the 
element  arrays  are: 


M|=  /  lA.iAjda  B|,y=  /  ^ikPjdrl  I 

jQe  Jrl  JQe 

Dy.,}  =  /  lAjdylAj 

JQe 


'j  B^x,ij  — 


L  = 


dr.  = 


^s,ij  - 


j  dxrpidxtl/jdQe  L‘yij  =  J  dyipjdyij/jdQe, 

f  dylpidxtl/jdQe  l-yx,o  =  [  dxtj/idyll/jdQe  ^  VrVsdPg, 

JQe  ’  JQe  J  fe 

f  v,v(dr,  f^=  f  v'XdCe  j 

Fe  J  Fe  J  Fg 

f  J/idxtpjdrl  Cljj=  f  >Pidytl/j^ 

jfI  JfI 


^jdrl 


/  lAjVrdrf, 

JFl 


(25) 


where  Qe  and  Fe  denote,  the  part  of  Q  and  F  associated  with  element  e.  Also,  ijji  and  V/  are  the  locally  defined  basis  functions 
from  which  the  global  basis  functions  ( Wi  and  t/)  are  formed.  For  quadrilateral  elements  with  spectral  order  N  as  used  in  this 
analysis,  will  be  discretized  into  {N+\f  points,  and  Vi  into  (N  +  1)  points,  thus,  ij  =  1,2,. .  .,(N+  1)^  and  r,  s  =  1,2,. .  .,N+  1. 
Now,  let: 


A  =  M-'  [Po^xBe  -  2UDx  -  2VDy] , 

B  =  [XxLx  +  XyLy  -fM  +  UVLxy  +  UVLyx  -  UVCn  +  UVCs  -  UVCe]  ,  (26) 

c  =  -XxM-^'aI 
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and  substitute  Eq.  (26)  into  Eq.  (22)  to  get  the  matrix  form  of  the  interior  problem: 

=  At/ +  Bf/ +  C(/)i.  (27) 

If  we  examine  the  boundary  auxiliary  variable  formulation  (13)  we  see  that  the  selection  of  appropriate  Q  values  for  the 
auxiliary  variables  has  not  yet  been  addressed.  It  should  be  noted  that  van  Joolen  et  al.  [31]  show  that  any  choice  of  Q  is 
guaranteed  to  reduce  spurious  reflection  as  the  order  of  the  NRBC  (J)  increases.  While  we  omit  the  details  here,  the  core 
of  this  argument  is  the  computation  of  a  so-called  reflection  coefficient  that  is  a  product  of  J  factors,  each  of  which  are  less 
than  one.  The  reflection  caused  by  the  artificial  boundary  must  decrease  as  the  order  of  the  NRBC  increases.  They  note,  “Of 
course,  a  good  choice  for  the  Q  would  lead  to  better  accuracy  with  a  lower  order  J,  but  even  if  the  ‘wrong’  C/s  are  taken. . .  one 
is  still  guaranteed  to  reduce  the  spurious  reflection  as  the  order  J  increases.”  [31,  p.  1045] 

Armed  with  this,  we  choose  convenient  values  for  our  C/s  that  cause  the  second  order  in  time  (o/)  terms  to  vanish.  In  the 
case  of  the  semi-infinite  channel,  on  the  eastern  boundary  this  value  is:  Q  =  Co  +  U.  A  physical  argument  for  this  choice  is  that 
since  the  predominant  speed  of  the  wave,  absent  the  advection  terms,  is  Co,  the  Q  terms  “correct”  the  boundary  formulation 
to  account  for  the  advection.  This  selection  for  the  C/s  then  makes: 

(x^  =  (X2  =  ■■■  =  a/_i  =  0, 

Ci=C2  =  ---  =  Cj-i=C, 

Now,  let: 

D=/2M'’-AyL'’  (28) 


and  substitute  (28)  into  (23)  to  get  the  following  form  of  the  problem: 

=  Dt]r  -  IXl’nr  +  +  XM'' ^2-. 

+  yD‘’<j>2  +  /IxMVs. 

=  002  +  +  2xM^04,  (29) 


i:P‘ii>j.2 + =  04>j_2 + 7d'’0j_i. 

If  we  now  collect  the  terms  on  the  left  and  right,  we  get  the  matrix  form  of  the  problem: 

Ed  =  ¥0  +  f],  (30) 

where: 


1  ^m'’ 

0 

0  \ 

/  yD" 

xM  0 

...  0  \ 

CD" 

0 

D 

...  0 

1  0 

0 

1  0 

0  0 

0  yDV 

(h  \ 
<l>2 

(  ] 
02 

/  D>7^  -  CD  V\ 

0 

0  = 

\  <l>]-\  / 

and  rj  = 

Vo/ 

It  should  be  noted  that  if  U  and  V  are  taken  to  be  zero,  that  this  entire  formulation  reduces  to  the  formulation  shown  by 
Lindquist  et  al.  [28]. 

5.  Time  integration 

The  formulation  outlined  in  (22)  and  (23)  constitute  a  system  of  coupled  ODEs  that  must  be  solved  to  yield  a  solution  for 
Since  the  goal  of  this  analysis  is  to  uncover  the  “true”  gains  made  by  high-order  boundary  treatment,  it  is  possible 
that  any  high-order  treatment  of  the  boundary  and  spatial  discretization  without  considering  a  high-order  treatment  of  the 
temporal  component  could  mask  gains  made  by  high-order  boundary  treatment.  Eor  this  purpose,  our  approach  uses  stan¬ 
dard  kth  order  Runge-Kutta  (RK)  methods  (up  to  order  10)  to  integrate  the  system  in  time. 

The  set-up  of  this  scheme  is  a  standard  one,  namely,  the  second  order  system  is  expanded  to  a  larger  system  of  first-order 
ODEs,  then  solved  appropriately  using  the  associated  RK  tableau.  Eor  most  cases  in  this  analysis  (unless  otherwise  stated) 
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time  integration  is  performed  using  a  4th  order  RK  scheme  using  a  time-step  chosen  to  ensure  a  Courant  number  of  0.25, 
where  the  Courant  number  is  defined: 

Courant  number  =  . 

Here,  ax  and  Ay  are  chosen  as  the  minimum  distance  between  any  two  points  in  the  x-  or  y-  directions  respectively.  Addi¬ 
tionally,  This  choice  is  made  since  the  interpolation  points  are  not  uniformly  distributed  when  using  spectral  elements. 


6.  Numerical  experiments 

Several  numerical  experiments  were  performed  to  solve  the  KGE  under  advection  with  and  without  dispersion.  In  order  to 
simplify  the  numerical  simulation  of  the  problem  at  hand,  the  KGE  (3)  is  converted  to  a  non-dimensional  form  as  described 
in  [32].  Using  typical  mesoscales  in  the  ocean,  the  length  scales  were  chosen  0(100  km),  vertical  depth  scales  0(100  m)  and 
the  dispersion  parameter /for  Coriolis  0(10“^  s“^).  In  terms  of  advective  velocities,  Majda  [32,  p.  61  ]  describes  typical  advec- 
tive  velocities  as  roughly  ^  that  of  the  medium  wave  speed  Co.  Here,  we  choose  to  challenge  the  boundary  further  by  choos¬ 
ing  faster  advective  velocities  at  ^  that  of  the  medium  wave  speed. 

Eor  experiments  that  follow,  the  reference  wave  speed  is  Co  =  1,  which  allows  the  initial  wave  to  propagate  through  the 
region  of  interest  in  a  reasonable  time  for  all  experiments.  Given  the  scale  choices  above  yields  a  dispersion  parameter/^  of 
0.1.  Eurther,  U  and  V  are  set  to  0.0  or  0.1  under  a  two-dimensional  Gaussian  initial  condition  to  test  the  formulation  in  a 
semi-infinite  and  infinite  channel  setting.  The  truncated  domain  Q  is  defined  on  x  =  [-2,2],  y  =  [-2,2]  and  the  initial  condi¬ 
tion  is  given  by: 

=  dtt](x,y,0)=0. 

In  order  to  see  the  effect  of  the  NRBC,  we  have  compared  our  solution  to  the  same  one  on  a  larger  domain  (x  =  [-2, 10]  for 
the  semi-infinite  channel,  and  x  =  [-6,6]  for  the  infinite  channel).  This  expanded  domain  reference  solution  replaces  the 
NRBC  with  no-fiux  boundary  conditions  as  well.  The  problem  is  then  solved  for  t  =  3  ensuring  that  the  disturbance  has  prop¬ 
agated  through  the  artificial  boundary,  yet  has  not  had  time  to  reach  the  no-fiux  boundary  and  pollute  the  reference  solution. 

To  quantify  the  errors  observed  between  the  reference  and  NRBC  solutions,  we  use  the  normalized  error  defined  as 
follows: 


error  ,2  =  , 


where  rj^  and  rj^  are  the  NRBC  and  reference  solutions. 

6.1.  Semi-infinite  channel 

Eor  the  first  experiment,  the  KGE  is  examined  on  a  semi-infinite  channel  with  the  NRBC  at  Xf  =  2  and  constant  advection 
with  U  =  0.1  and  \/=  0.  In  Eig.  5  we  plot  the  reference  solution  on  the  top  panel  and  the  solution  of  the  truncated  domain 


Expanded  Domain  Reference  Solution 


NRBC  Solution  (J=4) 


X 


Fig.  5.  Semi-infinite  channel  comparing  extended  domain  reference  solution  and  the  NRBC  solution.  4th  order  spectral  elements  using  NRBC  order  J  =  4 
with  advection  velocities  U  =  0.1,  V=  0.0  shown. 


Please  cite  this  article  in  press  as:  J.M.  Lindquist  et  al.,  Klein-Gordon  equation  with  advection  on  unbounded  domains  using  spectral  ele¬ 
ments  and  high-order  non-reflecting  boundary  conditions,  Appl.  Math.  Comput.  (2010),  doi:  10.1 01 6/j.amc.201 0.07.079 


ARTICLE  IN  PRESS 


10 


J.M.  Lindquist  et  al./ Applied  Mathematics  and  Computation  xxx  (2010)  xxx-xxx 


Fig.  6.  Differences  between  the  NRBC  solutions  [/=  1,5,10)  and  the  truncated  reference  solution  with  [7  =  0.1,  V=  0.0. 


using  the  G-N  NRBCs  on  the  bottom  panel.  The  solution  on  the  truncated  domain  uses  4th  order  spectral  elements  on  a 
24  X  24  element  grid,  discretizing  the  domain  into  9,409  points.  The  reference  solution  is  computed  on  the  extended  domain 
using  the  same  spectral  order  and  72  x  24  elements.  Qualitatively  speaking,  the  results  show  very  little  reflection  using  the 
combination  of  high-order  elements  and  NRBC.  The  differences  between  the  reference  and  NRBC  solutions  are  shown  in 
Fig.  6  for  NRBC  orders  of  1,  5  and  10. 

Quantitative  results  can  be  observed  in  Fig.  7  showing  the  error  on  ^2  as  a  function  of  spectral  and  NRBC  order 
(J=  1,. .  .,10,15,20).  The  number  of  elements  is  adjusted  for  each  spectral  order  to  maintain  an  equal  number  of  points 
(9,409)  that  the  domain  is  discretized  into.  It  is  clear  that  increasing  the  NRBC  order  yields  significant  gains  in  accuracy, 
but  by  7  =  5,  the  spatial  discretization  error  dominates  NRBC  error  in  the  low-order  element  (order  1  and  2)  cases.  However, 
by  increasing  the  spectral  order  of  the  elements,  this  spatial  discretization  error  decreases  and  allows  the  true  accuracy  of  the 
NRBC  to  be  realized.  The  Higdon  boundary  condition,  while  general  in  nature,  implicitly  assumes  that  by  the  time  the  wave 
pulse  gets  to  the  artificial  boundary,  it  is  traveling  primarily  as  a  plane  wave  normal  to  the  boundary.  The  previous  example, 
where  the  advective  velocity  was  in  the  same  direction  as  the  channel,  does  not  significantly  challenge  this  assumption.  In 
other  words,  to  really  test  the  value  of  the  boundary  condition,  one  must  try  advective  velocities  with  some  tangential  com¬ 
ponent  to  the  boundary.  Fig.  8(c)  and  (d)  show  the  same  plots  for  advection  velocities  in  other  directions,  one  being  the 
contrived  case  where  the  advection  is  perfectly  tangential  to  the  NRBC.  These  plots  correspond  with  the  contour  plots 
directly  below  that  show  the  comparative  solutions.  Examining  these  results,  it  is  clear  that  the  boundary  condition-even 
when  put  to  a  test  with  a  wave  pulse  containing  a  significant  tangential  component  to  the  boundary-performs  well.  It  is 
noted  that  the  order  of  the  error  suffers  slightly  under  diagonal  advection  when  compared  to  its  individual  axial  counter¬ 
parts.  This  may  be  due  to  the  additional  terms  activated  in  the  interior  and  boundary  formulations  (19)  and  (20)  when  U 
and  V  are  simultaneously  non-zero. 

6.2.  Infinite  channel 

For  the  next  set  of  experiments,  the  domain  is  an  infinite  channel  with  the  NRBCs  located  at  x  =  ±2.  The  set-up  is  similar  to 
that  of  the  semi-infinite  channel.  To  begin,  consider  the  Higdon  boundary  condition  for  the  western  boundary  given  by: 


Fig.  7.  Semi-infinite  channel  error  versus  NRBC  and  spectral  element  order.  Domain  is  discretized  into  9,409  points  for  all  spectral  element  orders  with 
advection  velocities  D=  0.1,  V=  0.0. 
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Expanded  Domain  Reference  Solution  Expanded  Domain  Reference  Solution 


(a)  r]{x,y,3):  U  =  0.0,  F  =  0.1 


(b)  v{x,y,3):  U  =  0.1,  V  =  0.1 


(c)  L^:  U  =  0.0,  V  =  0.1  (d)  L^:  U  =  0.1,  V  =  0.1 

Fig.  8.  Top  Plots:  Semi-infinite  channel,  4th  order  spectral  element  (J  =  4)  contour  plots  showing  rj{x,y,3)  on  extended  and  truncated  domains.  Domain  is 
discretized  into  9,409  points  for  all  spectral  element  orders  with  advective  velocities  specified.  Bottom  Plots:  Corresponding  error  versus  NRBC  and 
spectral  element  order. 


Expanded  Domain  Reference  Solution 


(a)  7?(x,  y,  3):  U  =  0.1,  V  =  0.1  (b)  L^:  U  =  0.1,  V  =  0.1 

Fig.  9.  Infinite  channel  comparing  extended  domain  reference  solution  and  the  NRBC  solution.  4th  order  spectral  elements  using  NRBC  order  J  =  4  with 
advection  velocities  Ll  =  0.1,  V=0A  shown. 


H, 


] 

n 

U=1 


f]  =  0  on  Fe. 


(31) 


Using  a  similar  strategy  for  introducing  a  set  of  auxiliary  variables  for  the  western  NRBC,  applying  them  to  the  KGE  equiv¬ 
alent,  then  converting  any  normal  derivatives  on  the  boundary  to  time  and  tangential  boundaries,  results  in  another,  very 
similar  formulation  that  is  directly  incorporated  into  the  weak  form  (18).  The  selection  of  parameters  Q  follows  the  “conve¬ 
nient”  choice  as  developed  for  the  eastern  boundary,  namely  to  remove  the  second  order  in  time  auxiliary  variable  term.  This 
choice  is  C,  =  Co  -  U. 


Please  cite  this  article  in  press  as:  J.M.  Lindquist  et  al.,  Klein-Gordon  equation  with  advection  on  unbounded  domains  using  spectral  ele¬ 
ments  and  high-order  non-reflecting  boundary  conditions,  Appl.  Math.  Comput.  (2010),  doi:  10.1 01 6/j.amc.201 0.07.079 


ARTICLE  IN  PRESS 


12 


J.M.  Lindquist  et  al./ Applied  Mathematics  and  Computation  xxx  (2010)  xxx-xxx 


Fig.  10.  Differences  between  the  NRBC  solutions  [/=  1,5,10)  and  the  truncated  reference  solution  with  [/=  V  =  0.1. 


Extended  Time  Solution 


-2-10  12 


Fig.  11.  Infinite  channel,  extended  time  solution  (t=  1000).  4th  order  spectral  elements  using  NRBC  order  J  =  4  with  advection  velocities  U  =  0A,  V=0.0 
shown. 


Similar  experiments  to  those  run  in  the  semi-infinite  channel  were  conducted  using  various  advective  velocities.  In  Fig.  9 
we  show  the  case  where  advection  is  towards  the  north-east  corner  of  the  domain.  Again,  qualitatively  and  quantitatively, 
the  NRBC  solution  performs  quite  well.  The  differences  between  the  reference  and  NRBC  solutions  are  shown  in  Fig.  10  for 
NRBC  orders  of  1,  5  and  10. 

6.3.  Long  term  results 

In  many  applications,  there  is  concern  for  the  stability  and  behavior  of  the  NRBC  scheme  in  the  long  term.  In  the  final 
infinite  channel  experiment,  we  have  run  the  KGE  with  7=10  using  4th  order  basis  functions  (9,409  global  points)  for 
t=  1000.  By  this  time  the  wave  front  has  exited  the  truncated  domain  and  the  reference  solution  is  expected  to  be  zero 
throughout  the  entire  domain.  The  differences  between  the  reference  (f/(x,y,1000)  =  0)  and  the  NRBC  solution  are  shown 
in  Fig.  11  are  of  order  10“"^.  This  figure  shows  that,  as  expected,  absent  sources  in  the  domain,  the  wave  front  exits  the  do¬ 
main  with  minimal  reflection  and  in  the  long  run,  tends  toward  a  steady  state  of  zero.  This  would  imply  long  term  stability  of 
the  method  when  used  in  this  context. 

6.4.  Effects  of  time  integration  technique 

At  the  outset  of  this  work,  it  was  believed  that  at  some  point  the  improvements  realized  by  improving  the  spatial  discret¬ 
ization  and  NRBC  would  eventually  be  limited  by  the  time  integration  scheme  [33].  To  this  end,  the  order  of  the  time  inte¬ 
gration  scheme  (RK2-RK10)  was  varied  to  examine  the  effects  of  time  integration  on  the  accuracy  of  the  solution.  As  has 
already  been  presented,  gains  made  by  increasing  the  order  of  the  NRBC  halt  for  lower  order  spectral  elements  after  7=5. 
Even  for  high-order  (order  8  and  16)  spectral  elements,  the  gains  made  by  increasing  the  order  of  the  NRBC  are  limited  at 
some  point  using  RK4.  In  [28]  the  authors  showed  that  high-order  time  integration  allowed  boundary  gains  to  improve  solu¬ 
tion  quality  for  the  KGE  under  zero  advection.  Here,  in  a  problem  incorporating  advection,  one  might  surmise  that  the 
changes  in  time  would  exaggerate  any  temporal  discretization  error. 

Eor  this  experiment,  consider  the  KGE  on  a  semi-infinite  domain  with  rj  =  0  on  Fw  To  ensure  that  any  boundary  or  time 
effects  are  not  masked  by  the  interior  discretization,  we  consider  8th  order  spectral  elements  discretized  into  9,409  global 
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Fig.  12.  Semi-infinite  channel  error  versus  RK  time  integration  order  using  various  order  NRBCs.  Domain  is  discretized  into  9,409  points  for  8th  order 
spectral  elements  with  advection  velocities  U  =  0.75,  17=  0.0  and  dispersion  parameter)^  =  1. 


points.  The  Gaussian  initial  condition  is  used  and  is  evaluated  until  t  =  3.  The  reference  solution  in  this  case  was  computed  as 
described  previously,  except  that  time  integration  was  performed  with  a  10th  order  Runge-Kutta  scheme  (RKIO)  using  a 
time-step  chosen  to  ensure  a  Courant  number  of  0.1.  Further,  to  exaggerate  the  issue  of  advection  and  dispersion,  we  con¬ 
sider  the  case  where  U  =  0.75  andf  =  1. 

As  can  be  observed  in  Fig.  12,  gains  made  by  improving  the  time  integration  matter  only  if  combined  with  high-order 
treatment  of  the  boundary.  Conversely-gains  using  high-order  treatment  of  the  boundary  can  only  be  realized  if  there  is  a 
sufficiently  high-order  treatment  of  the  time  integration.  While  this  analysis  conducted  time  integration  using  up  to  RKIO 
(18  stages),  RK4  or  RK5  appears  to  be  sufficient  in  exploiting  gains  when  using  high-order  boundary  treatment.  It  should 
be  noted  that  these  results  (error  on  the  order  of  10“^)  cannot  be  observed  unless  high-order  treatment  of  the  interior  also 
accompanies  the  high-order  treatment  of  the  boundary  and  time.  Several  experiments  were  conducted  which  varied  the  or¬ 
der  of  the  interior,  boundary  and  time  integration.  The  clear  result  was  that  without  high-order  treatment  of  all  components 
in  concert,  convergence  to  the  reference  solution  is  stalled. 

While  these  results  are  for  a  specific  problem  (KGE  with  advection),  we  believe  that  the  principle  of  a  balanced  approach 
to  all  components  (interior,  boundary  and  time)  is  a  sound,  extensible  procedure  for  any  problem.  If  high-order  treatment  of 
any  of  the  three  components  is  missing,  the  high-order  treatment  of  the  other  components  is  essentially  wasted. 

7.  Conclusions 

In  this  paper,  we  considered  a  reduced  form  of  the  shallow  water  equations  on  an  infinite  (channel)  domain.  Using  the 
Givoli-Neta  auxiliary  variable  formulation  of  the  Higdon  non-refiecting  boundary  conditions,  we  truncated  the  original  infi¬ 
nite  domain  and  developed  the  boundary  conditions  specific  to  the  problem  at  hand.  Using  a  high-order  approach  to  the  spa¬ 
tial  discretization  (spectral  element),  time  integration  (high-order  Runge-Kutta)  in  concert  with  the  high-order  boundary 
treatment,  we  showed  exponential  convergence  to  the  reference  solution  in  various  examples.  These  results  suggest  a  bal¬ 
anced  approach  to  dealing  with  truncation  errors  -  namely,  to  make  improvements  in  all  components  of  the  problem  to 
see  improved  accuracy. 

Future  work  includes  extending  these  results  to  the  quarter  and  open  plane  and  extending  all  of  this  work  to  the  full  first- 
order  SWE  system.  Challenges  in  extending  this  formulation  primarily  lie  in  developing  appropriate  corner  compatibility 
conditions  for  the  intersection  of  two  boundaries.  As  the  auxiliary  variable  formulation  is  itself  a  system  of  PDEs,  well-posed- 
ness  requires  appropriate  boundary  conditions  -  which  are  not  obvious  -  for  each  of  the  auxiliary  variables  when  employed 
in  a  “corner-type”  configuration.  The  authors  have  developed  stable,  low-order  boundary  conditions  for  the  auxiliary  vari¬ 
ables  for  the  KGE  with  advection,  albiet  with  less  impressive  results,  as  well  as  alternative  ways  of  handling  the  auxiliary 
variable  formulation  for  an  arbitrarily  shaped  convex  boundary  in  [34].  Work  is  ongoing  to  further  develop  the  auxiliary  var¬ 
iable  formulation  that  will  maintain  stable,  high-order  results  in  the  quarter  and  open  domain  experiments. 
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